Synchronization in disordered Josephson junction arrays: 
Small-world connections and the Kuramoto model 



B. R. Treefl and V. Saranathan 
Department of Physics and Astronomy 
Ohio Wesleyan University 

o ' 

<N ; D. Stroud 

Department of Physics 
, ^ \ The Ohio State University 

(Dated: February 2, 2008) 

Abstract 

We study synchronization in disordered arrays of Josephson junctions. In the first half of the paper, we consider the relation 
between the coupled resistively- and capacitively shunted junction (RCSJ) equations for such arrays and effective phase models 
of the Winfree type. We describe a multiple-time scale analysis of the RCSJ equations for a ladder array of junctions with 
non-negligible capacitance in which we arrive at a second order phase model that captures well the synchronization physics of 
the RCSJ equations for that geometry. In the second half of the paper, motivated by recent work on small world networks, 
we study the effect on synchronization of random, long-range connections between pairs of junctions. We consider the effects 
of such shortcuts on ladder arrays, finding that the shortcuts make it easier for the array of junctions in the nonzero voltage 
state to synchronize. In 2D arrays we find that the additional shortcut junctions are only marginally effective at inducing 
synchronization of the active junctions. The differences in the effects of shortcut junctions in ID and 2D can be partly 
understood in terms of an effective phase model. 
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I. INTRODUCTION 



The synchronization of coupled nonlinear oscillators has been a fertile area of research for decadesQ. In particular, 
phase models of the Winfree type 2] have been extensively studied. In one dimension, a generic version of this model 
for N oscillators is 

dft N 

^ = n j + J2^,kT(e k -e j ), (i) 
fc=i 

where dj is the phase of oscillator j, which can be envisioned as a point moving around the unit circle with angular 
velocity d6j/dt. In the absence of coupling, this overdamped oscillator has an angular velocity Vtj. T(9k — 0j) is the 
coupling function, and <Tj.k describes the range and nature {e.g. attractive or repulsive) of the coupling. The special 
case T(9k — 9j) = sin(9k — 9j), o~j.k = a/N (a = constant), corresponds to the uniform, sinusoidal coupling of each 
oscillator to the remaining N — 1 oscillators. This mean-field system is usually called the (globally-coupled) Kuramoto 
model (GKM). Kuramoto was the first to show that for this particular form of coupling and in the N — > oo limit, 
there is a continuous dynamical phase transition at a critical value of the coupling strength a c and that for a > a c 
both phase and frequency synchronization appear in the system (j, |4|. If o~j t k = a5j,k±i while the coupling function 
retains the form T(9j — Of.) = sin^, — 9j) we have the so-called locally-coupled Kuramoto model (LKM), in which 
each oscillator is coupled only to its nearest neighbors. Studies of synchronization in the LKM0, including extensions 
to more than one spatial dimension, have shown that a c grows without bound in the N — > oo limit |6(]. 

Several years ago, Watts and Strogatz introduced a simple model for tuning collections of coupled dynamical systems 
between the two extremes of random and regular networks 7] . In this model, connections between nodes in a regular 
array are randomly rewired with a probability p, such that p = means the network is regularly connected, while 
p = 1 results in a random connection of nodes. For a range of intermediate values of p between these two extremes, 
the network retains a property of regular networks (a large clustering coefficient) and also acquires a property of 
random networks (a short characteristic path length between nodes). Networks in this intermediate configuration 
are termed "small-world" networks. Many examples of such small worlds, both natural and human-made, have been 
discussed|8j- Not surprisingl y, t here has been much interest in the synchronization of dynamical systems connected 
in a small- world geometry pA Ho| . Generically, such studies have shown that the presence of small- world connections 
make it easier for a network to synchronize, an effect generally attributed to the reduced path length between the 
linked systems. This has also been found to be true for the special case in which the dynamics of each oscillator is 
described by a Kuramoto model [TTlll2] |. 

As an example of physically-controllable systems of nonlinear oscillators which can be studied both theoretically 
and experimentally, Josephson junction (JJ) arrays are almost without peer. Through modern fabrication techniques 
and careful experimental methods one can attain a high degree of control over the dynamics of a JJ array, and many 
detailed aspects of array behavior have been studied[13j. Among the many different geometries of J J arrays, ladder 
arrays (see Fig.^l deserve special attention. For example, they have been observed to support stable time-dependent, 
spatially-localized states known as discrete breathers 14] . In addition, the ladder geometry is more complex than 
that of better understood serial arrays but less so than fully two-dimensional (2D) arrays. In fact, a ladder can be 
considered as a special kind of 2D array, and so the study of ladders could throw some light on the behavior of such 2D 
arrays. Also, linearly-stable synchronization of the horizontal, or rung, junctions in a ladder (see Fig. ^) is observed 
in the absence of a load over a wide range of dc bias currents and junction parameters (such as junction capacitance), 
so that synchronization in this geometry appears to be robust [l5|. 

In the mid 1990's it was shown that a serial array of zero-capacitance, i.e. overdamped, junctions coupled to a 
load could be mapped onto the GKM 0,0. The load in this case was essential in providing an all-to-all coupling 
among the junctions. The result was based on an averaging process, in which (at least) two distinct time scales 
were identified: the "short" time scale set by the rapid voltage oscillations of the junctions (the array was current 
biased above its critical current) and "long" time scale over which the junctions synchronize their voltages. If the 
resistively-shunted junction (RSJ) equations describing the dynamics of the junctions are integrated over one cycle 
of the "short" time scale, what remains is the "slow" dynamics, describing the synchronization of the array. This 
mapping is useful because it allows knowledge about the GKM to be applied to understanding the dynamics of the 
serial JJ array. For example, the authors of Ref. ^(| were able, based on the GKM, to predict the level of critical 
current disorder the array could tolerate before frequency synchronization would be lost. Frequency synchronization, 
also described as entrainment, refers to the state of the array in which all junctions not in the zero- voltage state have 
equal (to within some numerical precision) time-averaged voltages: (h/2e)(d9j /dt)t, where 9j is the gauge-invariant 
phase difference across junction j. More recently, the "slow" synchronization dynamics of finite-capacitance serial 
arrays of JJ's has also been studiedfisl Il9j|. Perhaps surprisingly, however, no experimental work on JJ arrays has 
verified the accuracy of this GKM mapping. Instead, the first detailed experimental verification of Kuramoto's theory 
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was recently performed on systems of coupled electrochemical oscillators 20] . 

Recently, Daniels et al.[2l). with an eye toward a better understanding of synchronization in 2D JJ arrays, showed 
that a ladder array of overdamped junctions could be mapped onto the locally-coupled Kuramoto model (LKM). 
This work was based on an averaging process, as in Ref. 16], and was valid in the limits of weak critical current 
disorder (less than about 10%) and large dc bias currents, Ib, along the rung junctions (Ib/(Ic) ^ 3, where (I c ) is 
the arithmetic average of the critical currents of the rung junctions. The result demonstrated, for both open and 
periodic boundary conditions, that synchronization of the current-biased rung junctions in the ladder is well described 
byEq.m 

The goal of the present work is twofold. First, we will demonstrate that a ladder array of underdamped junctions 
can be mapped onto a second-order Winfree-type oscillator model of the form 

k=l 

where a is a constant related to the average capacitance of the rung junctions. This result is based on the resistively 
and capacitively-shunted junction (RCSJ) model and a multiple time scale analysis of the classical equations for the 
array. Secondly, we study the effects of small world (SW) connections on the synchronization of both overdamped and 
underdamped ladder arrays. We will demonstrate that SW connections make it easier for the ladder to synchronize, 
and that a Kuramoto or Winfree type model (Eqs. ^ an( i EJ, suitably generalized to include the new connections, 
accurately describes the synchronization of this ladder. 

This paper is organized as follows. In Sees. ITU and IIHI we discuss the multiple time-scale technique for deriving the 
coupled phase oscillator model for the underdamped ladder without SW connections. We compare the synchronization 
of this "averaged" model to the exact RCSJ behavior. We also analyze how the array's synchronization depends on the 
capacitance of the junctions. In Scc. lIVI we study the effects of SW connections, or shortcuts, on the synchronization 
of both overdamped and underdamped ladders. In our scenario, each SW connection is actually another Josephson 
junction. We generalize our phase-oscillator model to include the effects of shortcuts and relate our results to earlier 
work on Kuramoto- like models in the presence of shortcuts[lll[l3 • In Sec.Elwe study the effects of SW connections on 
synchronization in disordered 2D arrays. Here we find that the disordered 2D array, which does not fully synchronize 
in the pristine case (i.e. in the absence of shortcuts), is only weakly synchronized by the addition of shortcut junctions 
between superconducting islands in the array. In Sec. I VII we conclude and discuss possible avenues for future work. 



II. PHASE MODEL FOR UNDERDAMPED LADDER 



A. Background 

The ladder geometry is shown in Fig.^ which depicts an array with N = 8 plaquettes, periodic boundary conditions, 
and uniform dc bias currents, Ib, along the rung junctions. The gauge-invariant phase difference across rung junction 
j is 7j, while the phase difference across the off-rung junctions along the outer(inner) edge of plaquette j is ipi,j(ip2,j)- 
The critical current, resistance, and capacitance of rung junction j are denoted I C j, Rj, and Cj, respectively. For 
simplicity, we assume all off-rung junctions are identical, with critical current I co , resistance R a , and capacitance C Q . 
We also assume that the product of the junction critical current and resistance is the same for all junctions in the 
array 22], with a similar assumption about the ratio of each junction's critical current with its capacitance: 

IcjRj = IcoRo — (3) 



Icj_ _ Ico_ _ (Ic) , .% 

Cj ~ C ~ (C) ' [ 1 

where for any generic quantity X, the angular brackets with no subscript denote an arithmetic average over the set 
of rung junctions, (X) = (1/N) Y^Li X r 

For convenience, we work with dimensionlcss quantities. Our dimcnsionlcss time variable is 

t _ 2e(I c )t 
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FIG. 1: Ladder array with periodic boundary conditions and N — 8 plaquettes. A uniform, dc bias current Ib is inserted into 
and extracted from each rung as shown. The gauge-invariant phase difference across the rung junctions is denoted by jj where 
1 < j < A^, while the corresponding quantities for the off-rung junctions along the outer(inner) edge are r ^i,j{'^2,j)- The rung 
junctions are assumed to be disordered while the off-rung junctions are uniform. 



where t is the ordinary time. The dimensionless bias current is 

in * (6) 

while the dimensionless critical current of rung junction j is i c j — Icj / (-^c) ■ The McCurnber parameter in this case is 

_ 2e(/ c )(C) 

Note that (3 C is proportional to the mean capacitance of the rung junctions. An important dimensionless parameter 
is 

which will effectively tune the nearest-neighbor interaction strength in our phase model for the ladder. 

Conservation of charge applied to the superconducting islands on the outer and inner edge, respectively, of rung 
junction j yields the following equations in dimensionless variables: 

IB— icj smjj-i c j—^-i c:j fJ c ——^--asm'ipi j-a — ■r^--ap c — — r^+asm^ ,,-i+a p Vap c ri — = °> ( 9a ) 

dr dT l dr dr z dr dr z 



d^j d 2 li dip2 j n d 2 ip2 7 dip2 j-i n d 2 ip2 j-i 
-i B +i cj smj 3 +i cj — +i CJ {j c -^-asmip2,j-a—j^ — otp e rf 2 ' +asmip 2j j-i+ai — ^ Vap c — — = 0, (9b) 

where 1 < j < N. The result is a set of 2N equations in 3N unknowns: jj, ipi.j, and 4>2.j- We supplement Eq. [!]] 
by the constraint of fluxoid quantization in the absence of external or induced magnetic flux. For plaquette j this 
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constraint yields the relationship 



7j + ~>p2,j ~ 73+1 - ^1,3 = °- 



(10) 



Equations I§1 and llt)l can be solved numerically for the 3N phases jj, tpij and tbn psl] . 

We assign the rung junction critical currents in one of two ways, randomly or nonrandomly. We generate random 
critical currents according to a parabolic probability distribution function (pdf) of the form 



4A 3 



[A 2 - (ic - I) 2 ] , 



(11) 



where i c — I c / (I c ) represents a scaled critical current, and A determines the spread of the critical currents. Equationlll| 
results in critical currents in the range 1 — A<i c <l + A. Note that this choice for the pdf (also used in Ref. |lfj) 
avoids extreme critical currents (relative to a mean value of unity) that are occasionally generated by pdf 's with tails. 
The nonrandom method of assigning rung junction critical currents was based on the expression 
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[4? 2 -4(iV+l)j + (7V + l) 2 ] , l<j<N, 



(12) 



which results in the i C j values varying quadratically as a function of position along the ladder and falling within the 
range 1 — A < i c j < 1 + A. We usually use A = 0.05. 



B. Multiple time scale analysis 



Our goal in this subsection is to derive a Kuramoto-like model for the phase differences across the rung junctions, jj, 
starting with Eq.El We begin with two reasonable assumptions. First, we assume there is a simple phase relationship 
between the two off-rung junctions in the same plaquette: 



i>2,j = -4>i,j, 

the validity of which has been discussed in detail elsewhere pit |24| . As a result, Eg. 1101 reduces to 



1,3 



73 - 73 + 1 



(13) 



(14) 



which implies that Eg. I9al can be written as 



dr 2 



dr 2 



73 + 1 

dr 2 



- 2 



7j 



73-1 



dr 2 



dr 2 

ib — icj sin.7j +- a 



dr dr dr 



E 

<5=±1 



sin 



73+5 - 7j 



(15) 



Our second assumption is that we can neglect the discrete Laplacian terms in Eg 1151 namely V 2 [d^y j / 'dr) = 
djj + i/dr - 2d^.j/dr + d^fj-i/dr and V 2 (d 2 jj /dr 2 ) = d 2 ^ J+ i/dT 2 — 2d 2 jj/dr 2 + d 2 "/j-i/dr 2 . We find numerically, 



over a wide range of bias currents ig, McCumber parameters /3 C , and coupling strengths a that V {drfj/dr) and 
V 2 (c? 2 7j/(iT 2 ) oscillate with a time-averaged value of approximately zero. Since the multiple time scale method is 
similar to averaging over a fast time scale, it seems reasonable to drop these terms. In light of this assumption, Eo. 1151 
becomes 



d 2 ^ 

dr 2 



+ i 



dr/j_ 
C] dr 



IB 



H 3 



5=±l 



lj+S - 73 



We can use Eq. ^j]as the starting point for a multiple time scale analysis. Following Refs. 0] and \R 
Eg. 1161 bv is and define the following quantities: 



r = l B T 



(16) 



we divide 



(17a) 



$c = ibPc 



(17b) 
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e = l/i B - 



(17c) 



In terms of these scaled quantities, Eg. 1161 can be written as 

l = ioSc -^k + + ««* sin H - ™ E sin ( 7j+ V 7j ) ■ ( 18 ) 

Next, we introduce a scries of four (dimcnsionless) time scales, 

T n = e n T n = 0,1,2,3, (19) 

which are assumed to be independent of each other. Note that < e < 1 since e = I /is- We can think of each 
successive time scale, T n , as being "slower" than the scale before it. For example, T 2 describes a slower time scale 
than T\. The time derivatives in Eq. 1181 can be written in terms of the new time scales, since we can think of f as 
being a function of the four independent T n 's, f — t(Tq,Ti,T 2 ,T 3 ). Letting d n = d/dT n , the first and second time 
derivatives can be written as 

4= = d + eft + e 2 d 2 + e 3 8 3 (20) 



df 2 



8 2 + 2ed d 1 + e 2 (2d d 2 + d 2 ) + 2e 3 (d d 3 + d^) 



where in Eg. 1211 we have dropped terms of order e 4 and higher. 
Next, we expand the phase differences in an e expansion 



n=0 

Substituting this expansion into Eq. El an d collecting all terms of order e° results in the expression 

1 = i C j/3 c <9o7oj + *cjdo7oj, 

for which we find the solution 

T 

7o ,j = h 0j(Ti,T 2 ,T 3 ), 



(21) 



(22) 



(23) 



(24) 



where we have ignored a transient term of the form e~ Ta ^ c , and where (/)j(Ti, T2, T3) is assumed constant over the 
fastest time scale To- Note that the expression for 70 j consists of a rapid phase rotation described by To/i c j and 
slower-scale temporal variations, described by <f)j, on top of that overturning. In essence, the goal of this technique 
is to solve for the dynamical behavior of the slow phase variable, <f>j. The remaining details of the calculation can be 
found in the Appendix. We merely quote the resulting differential equation for the <pj here: 



d 2 4>j d(j) 



dr 2 



= Qj + Ki sin 



<t>j+6 - fa 



^E 



<5=±1 



5=±1 



i)j+S - fa 



bj+s ~ fa 
2 

<j+S ~ fa 



where Qj is given by the expression (letting Xj = i c j/iB for convenience) 



Qj = — 



and the three coupling strengths are 



(2/3 c 2 + : 



x) (3x 2 + 23/3 c 2 



16 (fi 2 c + x 



(25) 



(26) 



(27) 



6 



a x) (3/3 2 c - aj) 
icj 16 (/32+ x 2) 2 ' 



(28) 



M, 



(29) 



We emphasize that Eq. [53 is expressed in terms of the original, unsealed, time variable r and McCumber parameter 

We will generally consider bias current and junction capacitance values such that x 2 -C (3 2 . In this limit, Eas. 1271 
-|2!5]can be approximated as follows: 



,, CJ 



1 + 



1 



(30) 



a 



3^ 
16/3| 



O 



(31) 



Mi -> - — 




(32) 



For large bias currents, it is reasonable to truncate Eq. [25] at 0(1 /ig), which leaves 



/3 C 



dr 2 dr 



sin 



<5=±i 



6 J+<5 - <f>j 



(33) 



where all the cosine coupling terms and the third harmonic sine term have been dropped as a result of the truncation. 
In the absence of any coupling between neighboring rung junctions (a — 0) the solution to Eq. 1331 is 



i(a=0) 



= A + Be- T '^ + fli 



where A and B are arbitrary constants. Ignoring the transient exponential term, we see that 



,(«=o) 



/dr = Qj, so we 



can think of flj as the voltage across rung junction j in the uncoupled limit. Alternatively, flj can be viewed as the 
angular velocity of the strongly-driven rotator in the uncoupled limit. 

Equation 123 is our desired phase model for the rung junctions of the underdamped ladder. The result can be 
described as a locally-coupled Kuramoto model with a second-order time derivative (LKM2) and with junction coupling 
determined by a. In the context of systems of coupled rotators, the second derivative term is due to the non- 
negligible rotator inertia, whereas in the case of Josephson junctions the second derivative arises because of the 
junction capacitance. The globally- coupled version of the second-order Kuramoto model (GKM2) has been well 
studied; in this case the oscillator inertia leads to a first-order synchronization phase transition as well as to hysteresis 
between a weakly and a strongly coherent synchronized state|2a.l2q. 



III. COMPARISON OF LKM2 AND RCSJ MODELS 



We now compare the synchronization behavior of the RCSJ ladder array with the LKM2. We consider frequency 
and phase synchronization separately. For the rung junctions of the ladder, frequency synchronization occurs when 
the time average voltages, {vj) T — {d(j)j/dT} T are equal for all N junctions, within some specified precision. In the 
language of coupled rotators, this corresponds to phase points moving around the unit circle with the same average 
angular velocity. We quantify the degree of frequency synchronization via an "order parameter" 



/ = 1 



Sy(a) 
s v {0)' 



(34) 
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FIG. 2: (Color online) Frequency synchronization order parameter /, plotted versus nearest-neighbor coupling strength a for 
a ladder with TV = 10 plaquettes and bias current is = 5. Rung junction critical currents are assigned nonrandomly with 
A = 0.05. (a) j3 c = 1, (b) /3 C = 5, (c) /3 C — 10, (d) j3 c = 20. For each plot the phase differences and voltages are reset to zero 
with each new value of a. 



where s v (a) is the standard deviation of the N time-average voltages, 



u J/T . 



s (a) = \ 



2^=1 



N - 1 



(35) 



In general, this standard deviation will be a function of the coupling strength a, so s v (0) is a measure of the spread 
of the (vj) T values for N independent junctions. Frequency synchronization of all N junctions is signaled by / = 1, 
while / = means all N average voltages have their uncoupled values. 

Figure [21 compares the order parameter / for an array with N = 10 plaquettes, a bias current of is = 5, and 
nonrandomly-assigned critical currents with A = 0.05 for both the RCSJ model and the LKM2. For the RCSJ model, 
Eas. 151 and 1101 were solved numerically using a fourth-order Runge Kutta algorithm with a time step of At — 0.005 
and a total of 5 x 10 5 time steps. All time-average quantities were evaluated using the second half of the time interval. 
For the LKM2, the same numerical approach was applied to Eq. E31 Figure [3 shows some interesting behavior. First, 
in general, the LKM2 agrees well with the RCSJ model, especially in predicting a critical coupling strength, a c , at 
the onset of full frequency synchronization (/ = 1). Second, as (3 C is increased both models show evidence of a first 
order transition at a c (see Fig- Eld)) at which / jumps abruptly to a value of unity. In the vicinity of such an abrupt 
transition, the models differ the most, but even in Fig. Hfd), the RCSJ model and the LKM2 agree on the value of 
a c . The deviation between the models seen in Fig. E^d) near a ~ 0.4 could be due to a region of bistability near a c 
that becomes more prominent for increasing f3 c . 
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FIG. 3: (Color online) Frequency synchronization order parameter /, plotted versus nearest-neighbor coupling strength a for 
a ladder with N — 15 plaquettes, bias current is — 5, and /3 C = 20. Rung junction critical currents are assigned randomly with 
A = 0.025, and results are averaged over ten realizations of critical currents. Error bars represent the standard deviation of 
the mean value of / and for clarity only a few, representative error bars are shown. Numerical solution of the RCSJ model are 
denoted by filled circles, and the results from the LKM2 are denoted by crosses. 

Figure |3 shows the case where the critical currents are assigned randomly according to Eq. ^]with A = 0.025 
for N = 15, %b = 5, and (3 C = 20. The results for the frequency synchronization order parameter were obtained by 
averaging over ten different critical current realizations, and the error bars are the standard deviation of the mean 
value of / for each a. Note the excellent agreement between the RCSJ model and the LKM2. Also note that averaging 
over critical current realizations has a smoothing effect of / compared to, for example, Fig.|2Jd). 

Phase synchronization of the rung junctions is measured by the usual Kuramoto order parameter 



The results shown in Fig. ^represent the time-averaged modulus of r, (|r|) r , which approaches unity when the 
phase differences across the junctions are identical. Figure 0] compares the phase synchronization of the RCSJ model 
and the LKM2 for the same geometry as in Fig. [21 The agreement between the two models is excellent. Note the 
two types of behavior observable in the plots. For small coupling (a < 0.7), (|r|) T displays a complicated behavior 
due to finite-size effects, while for a > 0.7, (\r\) T exhibits a smooth rise toward a value of unity with increasing 
coupling. In fact, comparison of Figs. [21 and 0] shows that the value of a signaling the onset of the smooth increase 
in phase synchronization is approximately equal to a c , the value at which full frequency synchronization is obtained. 
Figure, ^d) also suggests that the finite-size fluctuations for small a are more pronounced at large j3 c (compare with 
Figs. Ufa), (b), (c)). Since the second-order Kuramoto model with global coupling (GKM2) has discontinuities in (\r\) T 
as a function of coupling strength for large arrays 25j|, and since we have mapped the RCSJ model to the LKM2, it 
would be interesting to look for evidence of a first-order transition in (|r|) T for large arrays. Such evidence is already 
visible, even for arrays as small as N — 10, in the frequency synchronization order parameter (see Fig. Hfd)). 

We have also studied the synchronization in our two models as a function of the dc bias current is for fixed coupling 
a, as shown in Fig.0 Such a graph is useful because experiments on periodic ladders would most likely be performed 
at fixed a (since that quantity is set by the fabrication of the rung and off-rung junctions), while the bias current 
could be easily varied. To obtain / experimentally, then, one needs to measure the time-average voltages across the 
rung junctions for each value of the bias current. Figure [5^ a) demonstrates that as the bias current is increased for 
fixed coupling strength, frequency synchronization is eventually lost. This is reasonable physically; as a rotator is 
driven harder a stronger coupling with its neighbors should be required to keep the rotators entrained. Figure [^b) 
plots s v (a, %b) versus is, showing that the spread in junction voltages scales linearly with the bias current over a wide 
range of currents. The behavior observed in both Figs. EJa) and (b) for bias currents of is ^ 10 is not surprising. 
When the system is far from frequency synchronization, the time-averaged voltages should be well approximated by 
their values in the absence of coupling, namely (vj) T « fij, where flj is given by Eq.[2<3 In the limit (i^/is) 2 < 1, 
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FIG. 4: (Color online) Phase synchronization order parameter (|r|) T , plotted versus nearest-neighbor coupling strength a for 
a ladder with N = 10 plaquettes, is = 5, and non-random critical currents with A = 0.05. (a) j3 c = 1, (b) f3 c = 5, (c) f3 c = 10, 
(d) /3 C = 20. For each plot phase differences and voltages across rung junctions are reset to zero with each new value of a. 



Eg. 1261 gives Qj ~ iB/icj- In this case, we can write 



Sti (a, is) 



IB 



\ 



/_L J_ 

N -1 



ClB, 



(37) 



where C is a constant independent of the bias current. Thus the linear scaling of s v with bias current is just due 
to the scaling of the time-averaged voltages across the rung junctions with %b- Equation 1371 is actually the standard 
deviation in the limit a — > 0, so for bias currents large enough that the junctions can be treated as approximately 
independent, we expect s v (a, iB)/s v (0, ib) — * 1, which in turn means / — > 0, as observed in Fig.|SJa). 

FigureElshows that ]im a —^o(vj) T = flj. To obtain this result, (v j) T (*s) across the j = 1 rung junction was calculated 
numerically for the RCSJ model for j3 c — 1 and a = 0.01, which is more than an order of magnitude smaller than 
a c . The results are shown as solid circles. The dotted line represents the analytic expression (vi) T = fii, where fii is 
given by Ea. 1261 and which results from our multiple time-scale analysis. The solid line is the large bias current limit 
of Eq.[2Hl namely w is/i c i. Note that the numerical results agree well with the Eq.JSHIfor a <C a c over the entire 
range of bias currents shown, and with the large bias current result for is Si 2.5. 

Of particular interest is how the array behaves near the frequency synchronization transition, a w a c . As shown 
in Fig. Efa) for an array with N = 10 plaquettes driven by a bias current is = 5, the order parameter develops a 
discontinuity at a c for (3 C 8. In addition, for sufficiently large (3 C the array also exhibits hysteretic behavior in 
/, as shown in Fig. The behavior depicted in Figs. [7| and is presumably due to bistability of the individual 
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FIG. 5: (Color online) (a) Frequency synchronization order parameter /, plotted versus dc bias current is for fixed coupling 
strength a — 0.25 for both the RCSJ model and the LKM2. N = 10, /3 C — 5, and nonrandom critical currents with A = 0.05. 
(b) Standard deviation of time-averaged rung junction voltages versus bias current for the same model parameters as in (a). 
The difference between the two models evident at large bias currents is probably due to phase slips in the off-rung junctions 
that violate Eg. 1131 



junctions[27j arising from their non-negligible capacitance. Figure[7fb) shows that for increasing j3 c the discontinuity 
in the order parameter at a c , A/, can be well fit by an exponential rise that asymptotically saturates to a value 
Mmax for C -> 00, i.e. 

i/ = {A/„., (38) 

For [] c < j3* the frequency synchronization transition is a smooth function of a as a decreases through a c from above. 
Thus, the junctions must be sufficiently underdamped for the discontinuous nature of the transition to be manifest. 
(But Fig.[7fa) gives at least the hint of a possible discontinuity in A/ for j3 c = 8 around a = 0.37 < a c .) The data in 
Fig. were obtained from a numerical solution of the RCSJ model, but we see qualitatively similar behavior from a 
numerical solution of Eq. 1331 namely saturation of A/ to a maximum value that is well fit by an exponential function. 

Lastly in this section, we address the issue of the linear stability of the frequency synchronized states (a > a c ) 
by calculating their Floquet exponents numerically for the RCSJ model as well as analytically based on the LKM2, 
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FIG. 6: Time- averaged voltage across j = 1 rung junction, plotted versus bias current is for TV = 10, a = 0.01 <C ot c , and 
j3 c = 1. The solid circles are a numerical result from the RCSJ model. They agree well with the analytical result for the 
weak-coupling limit, {vi) T = fii, where fl\ is given by Ea. 1261 which is represented by a dotted line in the graph. Deviations 
from the large bias current result, («i) T ~ is/ici (solid line), do not appear until is J$ 2.5. 



Eq. 123 The analytic technique used has been described in detail elsewhere 28j, so we shall merely quote the result 
for the real part of the Floquet exponents: 



Re(A m i c ) 



1 



l±ReJl-4p c (K + 3L) 



(39) 



where stable solutions correspond to exponents, A m , with a negative real part. One can think of the oj m as the normal 
mode frequencies of the ladder. We find that for a ladder with periodic boundary conditions and N plaquettes 



4 sin 2 (2^) 
l + 2sin 2 (^)' 



< m < N - 1. 



(40) 



To arrive at Eq. [3Hlwe have ignored the effects of disorder so that K and L are obtained from Eas. 1271 and 1281 with 
the substitution i c j — > 1 throughout. This should be reasonable for the levels of disorder we have considered (5%). 
Substituting the expressions for K and L into Eq. results in 



Re(A m i c ) 



1 

Wc 



1 ±Re 



\ 



1 - 2/3 c a { 1 



2/32 



(41) 



We are most interested in the Floquet exponent of minimum magnitude, Re(A mm t c ), which essentially gives the 
lifetime of the longest-lived perturbations to the synchronized state (see Fig- EJ - 

If the quantity inside the square root in Eq. ^2 is negative then Re(A m i c ) = — l/(2/3 c ). This is the value seen 
in Fig. EUa) for (3 C > f3 c = 1.56. For (3 C < f3 c and m = 1, the quantity inside the square root is positive and 

Re(A mm i c ) = (— l/2/3 c ) 1 — — 2f3 c auj\ , where we have used the fact that the quantity inside the braces in 

> 1. Physically, the crossover- type behavior evident in Fig.^Ja) is clue to 



Eq.^Jis essentially unity for is — 5 and (3 C 
the low frequency (long wavelength), m = 1, normal mode of the ladder changing from underdamped to overdamped 
in character as (3 C is decreased through (5 C = \j(2au>\) = 1.56 for = 10 and a = 1. Note from Fig. OJa) that 
the numerical result for the exponents based on the RCSJ model, with 5% disorder, agree quite well with Eq. ITTI 
Not surprisingly, the RCSJ model with no critical current disorder agrees very well with the analytic result since the 
disorder was ignored in order to obtain Eq. 1411 Figure OJb) shows how the minimum-magnitude Floquet exponent 
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FIG. 7: (a) (Color online) Frequency synchronization order parameter /, plotted versus coupling strength a for a ladder with 
N = 10, is = 5, and nonrandom critical currents with A = 0.05. Results based on the RCSJ model for ten different j3 c values 
are shown. In these simulations, a is initialized to a value greater than a c and then gradually decreased. With each new value 
of a the phase differences and voltages across the rung junctions are reset to zero. Note the appearance of a discontinuity 
in / as I3 C is increased. There is also some evidence of additional discontinuities in / for a < a c (see, for example, the data 
for /3 C = 8, 12). (b) Discontinuity in the frequency synchronization order parameter A/ at a — a c for different values of j3 c . 
This graph is produced from the data for / versus a seen in (a). The discontinuity is well fit by an exponential rise to a 
maximum, A/ = Af max [1 — exp — b(/3 c — /3*)], where we find parameter values of Afmax = 0.679 ± 0.020, b — 0.140 ± 0.020 
and P* = 7.750 ± 0.294 for this array. 

varies with coupling strength a for fixed (3 C . Now there is a crossover at a = a = l/(2j3 c ojf) = 1.56 for (3 C = 1. 
For a > a, Eg. fTTI gives Re(A mni i c ) = — l/(2/3 c ), independent of a. For a < a, however, we see that as a — » a c 
from above, the stability of the synchronized state decreases. In fact, one can show from Eq. ^2 that for a > a c and 
(3 c auj1 <C 1, the stability decreases linearly with a according to 

Re ( A min^J ~ g - , a -> aj 

independent of (3 C . Such linear behavior is evident in Fig.^b) for small a. 
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FIG. 8: Hysteresis in the frequency synchronization order parameter /, plotted as a function of the coupling strength a for a 
ladder with N = 10, ib = 5, /3 C = 30, and nonrandom critical currents with A = 0.05. To produce these plots, a is initialized 
to a value greater than q c and then gradually decreased until the system jumps discontinuously to the unsynchronized state; 
the coupling a is then increased until the system jumps back to the synchronized state. The final values of the phase differences 
and voltages across the junctions, as obtained from the previous value of a, are used as the initial values for each new value of 
a. (a) Results based on the RCSJ model, (b) Results based on the LKM2. Note that the models agree on the value of a c at 
which the jump occurs from the synchronized to the unsynchronized state for decreasing a. The solid lines are a guide to the 
eye. 

IV. "SMALL- WORLD" CONNECTIONS IN LADDER ARRAYS 

Many properties of small world networks have been studied in the last several years,including not only the effects of 
network topology but also the dynamics of the node elements comprising the network 8, 29]. Of particular interest has 
been the ability of oscillators to synchronize when configured in a small-world manner. Such synchronization studies 
can be broadly sorted into several categories. (1) Work on coupled lattice maps has demonstrated that synchronization 
is made easier by the presence of random, long-range connections |3Ct l31| . (2) Much attention has been given to the 
synchronization of continuous time dynamical systems, including the first order locally-coupled Kuramoto model 
(LKM), in the presence of small world connections [Til Il2l l32 ] | . For example, Hong and coworkers [Til have shown 
that the LKM, which does not exhibit a true dynamical phase transition in the thermodynamic limit (N — » 00) in 
the pristine case, does exhibit such a phase synchronization transition for even a small number of shortcuts. But the 
assertionj^ that any small world network can synchronize for a given coupling strength and large enough number 
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FIG. 9: (Color online) Real part of the minimum magnitude Floquet exponent for an array with N = 10 and is = 5. (a) 
Dependence of exponents on f3 c for fixed coupling strength, a = 1. Symbols are results of a numerical calculation based on the 
RCSJ model with either no disorder (open squares) or non-random critical currents and A = 0.05 (solid circles). The solid line 
is an analytic result (Eg. 1391 from the LKM2. (b) Dependence of exponents on coupling strength a for fixed (3 C = 1. 

of nodes, even when the pristine network would not synchronize under the same conditions, is not fully acce pted[3|. 
(3) More general studies of synchronization in small world and scale- free networks [ij have shown that the small 
world topology does not guarantee that a network can synchronize. In Ref. it was shown that one could calculate 
the average number of shortcuts per node, s sync , required for a given dynamical system to synchronize. This study 
found no clear relation between this synchronization threshold and the onset of the small world region, i.e. the value 
of s such that the average path length between all pairs of nodes in the array is less than some threshold value. 
Reference [T(il ] studied arrays with a power-law distribution of node connectivities (scale-free networks) and found 
that a broader distribution of connectivities makes a network less synchronizable even though the average path length 
is smaller. It was argued that this behavior was caused by an increased number of connections on the hubs of the 
scale-free network. Clearly it is dangerous to assume that merely reducing the average path length between nodes of 
an array will make such an array easier to synchronize. 

How do Josephson junction arrays fit into the above discussion? Specifically, if we have a disordered array biased 
such that some subset of the junctions are in the voltage state, i.e. undergoing limit cycle oscillations, will the 
addition of random, long-range connections between junctions aid the array in attaining frequency and/or phase 
synchronization? Our goal in this section of the paper is to address this question by using the mapping discussed 
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FIG. 10: A ladder array with periodic boundary conditions, N — 8 plaquettes, and one long-range connection. Rung junctions 
j = 1 and j = 4 are shown as connected by a pair of off-rung junctions, which have phase differences of ^1:1,4 an d i>2;i,4- The 
additional off-rung junctions are uniform and identical to the off-rung junctions in the pristine ladder. Uniform dc bias currents 
are applied along the rungs as in Fig. 

in Sees. ITT1 and IIHI between the RCSJ model for the underdamped ladder array and the second-order, locally-coupled 
Kuramoto model (LKM2). Based on the results of Ref. |2l|, we also know that the RSJ model for an overdamped 
ladder can be mapped onto a first-order, locally-coupled Kuramoto model (LKM). Because of this mapping, the ladder 
array falls into category (2) of the previous paragraph. In other words, we should expect the existence of shortcuts 
to drastically improve the ability of ladder arrays to synchronize. 

We add connections between pairs of rung junctions that will result in interactions that are longer than nearest 
neighbor in range. We do so by adding two, nondisordered, off-rung junctions for each such connection. For example, 
Fig. 1101 shows a connection between rung junctions j = 1 and j — 4. This is generated by the addition of the two 
off-rung junctions labeled tpi-,iA arL d ^2:1,4, where the last two indices in each set of subscripts denote the two rung 
junctions connected. The new off-rung junctions are assumed to be identical to the original off-rung junctions in the 
array, with critical current I co , resistance i? D , and capacitance C for the underdamped case. Physically, we should 
expect that the new connection will provide a sinusoidal phase coupling between rung junctions j = 1 and j ' = 4 with 
a strength tuned by the parameter a = I co /(I c ), where (J c ) is the arithmetic average of the rung junction critical 
currents. We assign long range connections between pairs of rung junctions randomly with a probability p distributed 
uniformly between zero and one, and we do not allow multiple connections between the same pair of junctions. For 
the pristine ladder, with p = 0, each rung has only nearest-neighbor connections, while p = 1 corresponds to a regular 
network of globally coupled rung junctions, i.e. each rung junction is coupled to every other rung junction in the 
ladder. 

In Fig. ^2 we plot the two standard quantities used to characterize the topology of the network: the average path 
length I and the cluster coefficient C, calculated numerically for a network with N = 50 nodes (i.e. rung junctions). 
The average path length / is defined as the minimum distance between each pair of nodes averaged over all such 
pairs, while the cluster coefficient C is the average fraction of nodes neighboring each node that are also neighbors 
themselves. These quantities are plotted as a function of the product pN. For pN ps 1, the average path length is 
already substantially reduced from its value in the pristine limit, p = 0. It is this reduced average distance between 
pairs of nodes that is one of the hallmarks of a small-world network. Because our ladder array, in the pristine limit, 
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FIG. 11: Scaled average path length l/N (solid circles and left vertical axis) and cluster coefficient C (hollow circles and right 
vertical axis), plotted versus the product pN, where p is the probability of making a shortcut connection between pairs of nodes 
and N is the total number of nodes in the network. This graph corresponds to N = 50. The solid line is C = 2/JV, which 
is the expected cluster coefficient for a random network of size N = 50 in which each node has two neighbors in the pristine 
limit. Note that C — > as pN —* 0, since each node's neighbors are not connected to each other in that limit. Also note that 
C — * 1 as p — * 1, which is the limit of a fully-connected, regular network in which each node's neighbors are all connected to 
one another. 



allows only nearest-neighbor coupling, the cluster coefficient vanishes as p — > 0. As a result, our ladder geometry 
does not conform to the most commonly-accepted definition of a small world, in which both reduced path lengths 
(compared to the pristine limit) and high cluster coefficients (compared to that of a random network) coexist (see 
Fig. 2 in Ref. @). Nevertheless, in our simulations we routinely choose value for the parameter p such that pN places 
us in the region of reduced path lengths. So we are considering ladder arrays in which the average distance between 
rung junctions is reduced by shortcuts by a factor of five to ten. 

Next, we argue that the RCSJ equations for the underdamped junctions in the ladder array can be mapped onto 
a straightforward variation of Eq. 1331 in which the sinusoidal coupling term for rung junction j also includes the 
longer-range couplings due to the added shortcuts. Imagine a ladder with a shortcut between junctions j and I, where 
I ^ j,j ± 1. Conservation of charge applied to the two superconducting islands that comprise rung junction j will 
lead to equations very similar to Eq. |5J For example, the analog to Ea. l9al will be 

i B - i C j sin.7j - icj-j^ ~ Pcicj^f - asin^ij - a ~~^ L ~ Pc a , ^ + (42) 



a sin»i + a h p c a — 5 h > asinfei + a — ; \- p c a — — 5— 

dr dT A < 1 dr dr A 



0. 



with an analogous equation corresponding to the inner superconducting island that can be generalized from Eq. 19 bl 
The sum over the index I accounts for all junctions connected to junction j via an added shortcut. Fluxoid quantization 
still holds, which means that we can augment Eq. 1101 with 

Jj + 4>2;jl ~ li ~ ipi;jl = 0. (43) 
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We also assume the analog of Eg. 1131 holds: 



1p2;jl = -ipi;ji- (44) 
Equations 1431 and 1441 allow us to write the analog to Eq. ^]for the case of shortcut junctions: 

fctfl = (45) 

Equation 1431 in light of Eq. 021 can be written as 

i B - u BtaTi - ^ - teA + «E sin f +«E sin f ^) + ( 4 6) 



2 Urj 2 ^J + 2^U rfr / 2 Y ^ dr 2 dr 2 

where the sums £; are over all rung junctions connected to j via an added shortcut. As we did with the pristine ladder, 
we will drop the two discrete Laplacians, since they have a very small time average compared to the terms i C jd'jj/dT + 
i C jf3 c d 2 jj /dr 2 . The same is also true, however, of the terms a/2 ^li^lj / dr—d^i/ dr) and a/2 J2i(d 2 Jj/dT 2 —d 2 ^i/dT 2 ), 
as direct numerical solution of the full RCSJ equations in the presence of shortcuts demonstrates (see Fig. ll2[l . So we 
shall drop these terms as well. Then Eq. ^\ becomes 

i B - i cj sin 7j ~ i cj ^ - c i cj — §L + | E sin i 1 -^ 21 ) ' ( 4? ) 

where the sum is over all junctions in Aj, which is the set of all junctions connected to junction j. Based on our work 
in Sec.|ni we can predict that a multiple time scale analysis of Eq. 07| results in a phase model of the form 

^ + f = O j + f Esi „(^4), (4S> 

fceA 3 v 7 

where fij is give by Eq. 1261 A similar analysis for the overdamped ladder leads to the result 

keAj v 7 

where the time-averaged voltage across each overdamped rung junction in the uncoupled limit is 



Of) = 1 /fg) 2 -l. (50) 



Figure ITSI demonstrates that the frequency synchronization order parameter /, calculated from Eq . I49I for over- 
damped arrays with A = 30 and N = 50 and in the presence of shortcuts, agrees well with the results of the RSJ 
model. In addition to the pristine array with p — 0, we considered arrays with p = 0.05 and p = 0.10 in which we 
averaged over 10 realizations of shortcuts. The agreement between the two models is excellent, as seen in the figure. 
It is also clear from the figure that shortcuts do indeed help frequency synchronization in that a smaller coupling 
strength a is required to reach / = 1 in the presence of shortcuts. In fact, the value of a c required to reach / = 1 
is growing with increasing N in the cases of p = 0; for example, we find the a c = 3.36 for N — 100 (compare with 
Fig. |13|l. The same is clearly not true, however, for arrays with p = 0.05 and p = 0.1. Recently, Hong, Choi, and 
Kim|ll| have demonstrated, using a finite-size scaling analysis applied to the LKM, that the phase synchronization 
order parameter, (|r|) T , in the presence of shortcuts (0 < p < 1) has a mean-field synchronization phase transition as 
in the GKM (i.e. LKM with p = 1). Based on the agreement between the two models shown in Fig. ^| we expect 
such an analysis to apply to the RSJ equations for the ladder as well. The finite-size scaling behavior of the frequency 
synchronization of the LKM has not been studied, so the nature of that transition is not well known. Figure ITU 
demonstrates that underdamped ladders also synchronize more easily with shortcuts and that Eq. 1481 agrees well with 
the RCSJ model. 

Although the addition of shortcuts makes it easier for the array to synchronize, we should also consider the effects 
of such random connections on the stability of the synchronized state. The Floquet exponents for the synchronized 
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FIG. 12: Time dependence of several combinations of voltages or voltage derivatives for a ladder with N = 10, is = 5, and 
nonrandom critical currents with A = 0.05. The array was chosen to have three shortcuts between the following pairs of rung 
junctions: (1,3), (2,7), (4,6). In both plots the following quantities are compared (note that Vi = dfi/dr): i c iVi + i c i/3 c dvi/dr 
(dot-dashed line), (a/2)(wi — 113) (solid line), and (a/2)/3 c (dvi/dT — dv^/dr) (dashed line). The time average value of the latter 
two quantities is negligible, (a) a = 1, f3 c = 10. (b) a = 1, j3 = 1. 



state allow us to quantify this stability. Using a general technique discussed in Ref. |35| , we can calculate the Floquet 
exponents X rn for the LKM based on the expression 

\ m t c = aEg, (51) 

where are the eigenvalues of G, the matrix of coupling coefficients for the array. A specific element, Gij, of this 
matrix is unity if there is a connection between rung junctions i and j. The diagonal terms, Gu, is merely the negative 
of the number of junctions connected to junction i. This gives the matrix the property G\j — 0. In the case of 
the pristine ladder, the eigenvalues of G can be calculated analytically, which yields Floquet exponents of the form 

A m (^ c = -4asin 2 (^). (52) 

This result is plotted in Fig. ^j] as the solid line for an overdamped array with N ~ 100; note that the solid line 
is the p = Floquet exponent of minimum, nonzero magnitude. Since the are purely geometry dependent, i.e. 
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FIG. 13: (Color online) Frequency synchronization order parameter /, plotted versus coupling strength a for overdamped 
arrays with bias current i s = 5 and non-random critical currents with A = 0.05. Solid symbols represent numerical results 
based on the LKM model, Ea. 1491 while hollow symbols are from the RSJ model. The data for nonzero p represent an average 
over ten realizations of randomly-assigned shortcuts, (a) N — 30, (b) N = 50. Note that shortcuts improve the frequency 
synchronization behavior in that the critical coupling needed for / = 1 is clearly reduced as p is increased from zero. In fact 
a c is clearly growing with increasing N in the case of p = 0, while the same is not true of the array with shortcuts. 



do not depend on the coupling strength, we expect the exponents to grow linearly with a, based on Eq. 1521 To 

include the effects of shortcuts, we found the eigenvalues numerically for a particular realization of shortcuts 

(for a given value of p), and then we averaged over 100 realizations of shortcuts for each value of p. The exponents 

of minimum magnitude for the overdamped array with p = 0.01,0.05,0.10 are also shown in Fig. 1151 (note the 

logarithmic scale on both axes). Clearly, shortcuts greatly improve the stability of the synchronized state. Specifically 

^(p— o.i) ,^(p—0) _ 2030 a three-order of magnitude enhancement, 
mm ' mm ' ° 
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FIG. 14: (Color online) Frequency synchronization order parameter /, plotted versus coupling strength a for an underdamped 
ladder array with TV = 30, bias current is = 5 and non-random critical currents with A = 0.05. Solid symbols represent 
numerical results based on the LKM2 model, Eq. 1481 while hollow symbols are from the RCSJ model. The data for nonzero 
p represent an average over ten realizations of randomly-assigned shortcuts, (a) (3 C — 1, (b) /3 C == 30. Agreement between the 
two models is excellent, except for large /3 C and p = near a c , which is not unexpected based on Fig. 

V. "SMALL- WORLD" CONNECTIONS IN TWO-DIMENSIONAL ARRAYS 

In this section we present some preliminary results on synchronization in disordered two-dimensional (2D) arrays 
in the presence of shortcuts. Geometrically, we can think of a pristine 2D array as a set of M columns, or "ladders" , 
each with N plaquettes, grafted together (see Fig. which depicts M = 2, N = 5). It is well known that in 
such a geometry, phase locking of all the horizontal junctions can occur in a horizontally-biased uniform array (i.e. 
no critical-current disorder) but that a high-degree of neutral stability is exhibited|36j|. More precisely, in an array 
with M columns, there will be a zero-valued Floquet exponent with multiplicity M. It was shown in Ref. 37] that 
underdamped arrays in an external magnetic field perpendicular to the plane of the array could lift this degeneracy via 
the coupling of the junctions to the external magnetic field in the gauge-invariant phase difference. In the presence of 
critical current disorder, however, numerical simulations of the RSJ equations revealed that frequency synchronization 
was no longer possible. In fact, each column or ladder in the array would individually synchronize but that sufficient 
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FIG. 15: Floquet exponent of minimum magnitude A mm , plotted versus coupling strength a for the LKM with N = 100. 
Exponents are calculated based on a technique described in Ref. |3S||. For p = 0, the exponents can be calculated analytically, 
with the result Re(A m : n t c ) = —4a sin 2 (tt/N). For each nonzero p, the results are averaged over 100 realizations of shortcuts. 
Note the logarithmic scale on both axes. 



"inter-ladder" coupling was not present to entrain the entire set of horizontal junctions 38, 39] . This behavior is shown 
in Fig. El by means of a so-called cluster diagram. 

To include shortcut connections in the 2D array, we followed a procedure similar to that described in the previous 
section. For clarity, let each horizontal junction in the array be described by a pair of coordinates where i denotes 
the row and j the column in which the junction is positioned. To establish a connection between two particular 
horizontal junctions, say junctions i,j and k,l, that are not already nearest neighbors, we add two new junctions, one 
connecting the left superconducting island of junction i, j to the left island of junction k, I and the second junction is 
added between the islands on the right sides of i, j and k, I. We also allow the shortcut junctions to be critical-current 
disordered. Contrary to the case of individual ladder arrays, the effects of shortcuts on synchronization in 2D arrays is 
not so easily characterized. Figure ITS1 shows the scaled standard deviation of the time-averaged voltages, s v (a)/s v (0), 
for the ten horizontal junctions of an overdamped array with M — 2 and N — 5 and for p — 0.25. Note that the 
ratio does not approach unity for nonzero p as a — > because of the presence of the disordered shortcut junctions. 
For reference, s v (a) / s v (0) for the pristine array is also shown (solid circles). In this case, entrainment is frustrated in 
that the ratio settles into a clearly nonzero value as the coupling strength is increased. The hollow circles and squares 
in Fig.^|are the values of s v (a) / s v (0) for two different realizations of shortcuts at p = 0.25. For case 2 in the Figure 
(hollow squares), the shortcuts have only slightly improved the level of synchronization compared to the pristine case, 
as s v (a)/ s v (0) is only reduced by a factor of about 0.35 compared to its p — value. Case 1 (hollow circles) is more 
interesting in that s v (a)/s v (Q) is reduced, on average, by an order of magnitude compared to the pristine case by the 
particular realization of shortcut junctions present. (Note the logarithmic scale on the vertical axis and the topmost 
arrow on the right axis, which denotes the average value of s v (a) / s v (Q) for 4 < a < 10.) The noise evident in the 
results for case 1 is probably a finite-size effect, but studies of larger arrays are necessary to be sure. 

Although array synchronization has clearly not occurred in the second realization of shortcuts in Fig. 1181 the 
reduced value of s v (a)/s v (0) for case 1 does not automatically imply entrainment has occurred in that case. Included 
in Fig. 1181 are the values of s v (a)/s v (0) for the junctions in each column separately (hollow triangles). The low average 
value of these quantities (see the two lower arrows along the right axis) shows that the junctions in a given column 
are much more strongly entrained to each other than to junctions in the neighboring column. Thus, the hollow circles 
in Fig. 1181 correspond, at best, to weak intercolumn synchronization. Figures Il9f a) and (b) are cluster diagrams for 
cases 1 and 2, respectively, in which the vertical axes of the two plots have the same scale. Simple visual inspection 
of the plots suggest that the array is weakly frequency synchronized in case 1 for a > 2 but not in case 2. In fact, we 
have considered ten different realizations of shortcuts at p — 0.25 and find this weak synchronization behavior only 
for one realization (i.e. case 1). The nine remaining realizations resulted in an array that was clearly not entrained, 



22 



1. 1 

-X— 



1,2 

-X— 



>: 



>: 



+y 



>: 



+x 



1. 1 



1,2 



FIG. 16: A two dimensional array of junctions with M — 2 columns, or "ladders", each with N — 5 plaquettes. A dc bias 
current Ib is injected along the left side and extracted along the ride side of the array. We assume periodic boundary conditions 
in the y direction. A junction with a label such as 1, 2 means the junction is in row 1 and column 2. 

as in case 2. Based on these results, we thus conclude that shortcuts in 2D arrays, biased in the standard way shown 
in Fig. El do not significantly enhance the array's ability to synchronize. We see similar effects for p = 0.5. 

We have also considered the effects of uniform shortcuts in the 2D array, where each additional shortcut junction 
is identical to the uniform vertical junctions in the pristine array. As in the case of disordered shortcuts, when we 
consider ten different realizations at p — 0.25 we find that in some cases (roughly 40% of the realizations) the array 
weakly frequency synchronizes and in the remaining cases it clearly does not. Results representative of these two 
outcomes are shown in Fig. [2U] They show an interesting distinction between this case and the case of disordered 
shortcuts discussed previously: for sufficiently large coupling a, all average voltages across the horizontal junctions 
now go to zero. This behavior is due to the fact that a — I C o/(Ic) is the ratio of the critical current of the vertical 
junctions, now including shortcut junctions, to the average critical current of the horizontal junctions. As this ratio 
increases, for a given number and configuration of shortcuts, a value of a is eventually reached for which all the bias 
current is able to traverse the circuit without exceeding any particular junction's critical current. In other words, the 
array is biased below its effective critical current in the presence of shortcuts and thus there is a zero average voltage 
across the array. 

The limiting case of p — 1 means that each horizontal junction is connected to all the remaining horizontal junctions. 
For this case and with uniform shortcut junctions, we find that the array behaves similarly to Fig. I20f a): there is a 
range of coupling strengths a for which there is weak entrainment, but for large a the array is in the zero- voltage 
state (see Fig.l2T]). 



VI. CONCLUSIONS 



In this paper we have obtained two main sets of results. First, using a multiple time scale method, we have mapped 
the exact RCSJ equations for an underdamped ladder with periodic boundary conditions to a second order, locally- 
coupled Kuramoto model (LKM2). Secondly, we have studied the effects of small world connections on the ability of 
both ladder and 2D arrays to synchronize. The mapping to the LKM2 is itself useful for two main reasons. First, the 
synchronization behavior of the Kuramoto model and its variations has been well studied in its own right and could 
thus shed some light on behavior of actual JJ arrays. Secondly, the LKM2 is solved more quickly on the computer 
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FIG. 17: A cluster diagram of the time-averaged voltages, (vij) T , across the horizontal junctions, plotted versus coupling 
strength a for an array with M — 2 columns and N = 5 plaquettes per column. A bias current is = 5 is applied and extracted 
along each row, and the horizontal junctions are assigned critical currents randomly according to Ea. 1111 with A = 0.05. We 
assume periodic boundary conditions in the vertical direction. The array is pristine (p = 0). The symbols with dots correspond 
to voltages across the junctions in column 2. The diagram demonstrates that the two columns are each frequency synchronized 
but that the synchronized voltages for each column are different. The numbers in the legend denote the coordinates of each 
horizontal junction. For example the coordinates 1, 2 denote the horizontal junction in the first row (starting from the bottom 
row) and second column (starting from the left column). 



and is easier to understand intuitively than the RCSJ equations. Future work in this area could include using the 
results of this mapping for the underdamped ladder (as well as the first-order LKM for the overdamped ladder) to 
arrive at a phase model for the 2D array, as has been suggested in Ref. [40|. Such a phase model for 2D arrays may 
shed light on why it is difficult for 2D arrays to synchronize. 

We have also shown that small-world connections enhance a ladder array's ability to synchronize. This result is not 
surprising in light of our mapping to the LKM and earlier studies in which the LKM was found to exhibit a mean- 
field like phase synchronization transition in the presence of shortcuts But we find that SW connections only 
marginally increase the ability of a 2D array to synchronize. Specifically, for several representative 2D small-world 
networks, we found only some small fraction of these networks had slightly improved frequency synchronization of the 
horizontal junctions. In the pristine 2D array (p = 0) it is well known that no such synchronization, weak or otherwise, 
is observed over the entire array- hence our characterization that shortcuts are only marginally effective at producing 
a synchronized array. This conclusion holds whether the additional shortcut junctions are disordered or uniform. 
Future work in this area could include looking at a broader range of p values as well as looking at larger 2D arrays. It 
is tempting, but an oversimplification, to think of 2D arrays as mere assemblages of ladder arrays. One source of this 
temptation is the intriguing fact that the pristine 2D array of disordered junctions will form synchronized clusters 
consisting of individual ladders. Another is our result that shortcuts can augment synchronization in individual 
ladders but not in 2D arrays. If one can produce a mapping, even approximate, of the RCSJ equations for the 2D 
array to a phase model of the Winfree type, this could be very helpful in understanding the rich and perplexing 
dynamical behavior of 2D arrays. 
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FIG. 18: (Color online) Scaled standard deviation s v (a) / s v (0) of time-averaged voltages across the horizontal junctions, 
plotted versus coupling strength a for an array with M = 2 columns and N — 5 plaquettes per column. The bias current is 
is = 5, and the horizontal junctions are assigned critical currents randomly according to Eg. 1111 with A = 0.05. We assume 
periodic boundary conditions in the vertical direction. The solid circles are for the pristine array, p = 0. The hollow circles and 
squares represent two different realizations of shortcuts at p — 0.25 in which the added shortcut junctions are critical-current 
disordered. The hollow up(down) triangles denote the value of s v (a)/s v (0) for only the junctions in column one(two). The 
arrows pointing to the right axis denote the average values of s„(a)/s„(0) over the interval 4 < a < 10 for the hollow circles 
and both the up and down triangles. 



APPENDIX: MULTIPLE TIME SCALE ANALYSIS 



We substitute Eq.|22into Eq.^J To help organize the terms according to the order in e, we also write the nonlinear 
terms as follows (as in Ref. (I9j1: 

oc 

sin7j = e n Sn,j (A. la) 

n=0 
oo 

= j2 enR ^+s> ( A - lb ) 

where Soj = sin7o, 3 -, Sij = 7ijCOS7oj, S 2 .j = 72jCOS7 0j - §7? tj sin 70, j , Ro,j+s = sin [(70,3+a - 7o,j)/2], 
Rij+6 = i cos [(7o,j-H5 -7o,j)/2] (7ij+«5 - 7ij), and R 2 j+s = § cos [(-y ,j+6 - 7o,j)/2] (72,3+* ~ 72,j) - 
i sin [(7o,j+5 — 7o.j)/2] (~fi,j+s — 7i,j) 2 - So Eq.^|can then be written as 

00 

1 = icjfic e " i d o + 2ed ° d i + e " ( 25 o9 2 + d\) + 2e 3 (9 a 3 + frfc)] 7 «,j + 

00 oc 00 

£ e " + e9 i + ^ + ^3] 7 n,3 + « CJ J] e " 5 n,3 " e " £ E eW ^+«- ( A - 2 ) 
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FIG. 19: Cluster diagrams for the time-averaged voltages, {vij) T , across the horizontal junctions, plotted versus coupling 
strength a for an array with M — 2 columns and N = 5 plaquettes per column. The bias current is is =5, applied 
and extracted along each row, and the horizontal junctions are assigned critical currents randomly according to Eq. 1111 with 
A = 0.05. We assume periodic boundary conditions in the vertical direction. The legend provides the coordinates for each 
horizontal junction as in Fig. 1171 (a) p — 0.25, first shortcut realization, (b) p — 0.25, second shortcut realization. 



Extracting all terms of O(e ) yields 



1 = i c i/3 c <9o7o,j + icjd "fo,r 



(A.3) 



The solution to the homogeneous version of Eq. lXTSl is jqj = A + Be T o/0<= : where A and B are constants with respect 
to the Tq time scale. The exponential term is dropped because it represents transient behavior. We take a particular 

solution to the inhomogeneous equation of the form = C^Tq, where cj ^ is independent of Tq. Substitution of 
into Eg. I A. 31 yields cj ' = l/i c j- So the solution to Eg. IA. 31 can be written as 



7o j = h 4>j(Ti,T 2 ,T 3 



(A.4) 



where one can think of <j)j as describing the slow phase dynamics of 70, 
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FIG. 20: Cluster diagrams for the time averaged-voltages, {vij) T , across the horizontal junctions, plotted versus coupling 
strength a for an array with M = 2 and N = 5. The bias current is is = 5, and the horizontal junctions are assigned critical 
currents randomly according to Eg. 1111 with A = 0.05. We assume periodic boundary conditions in the vertical direction. 
All shortcut junctions are identical to the uniform vertical junctions in the pristine array, i.e. the shortcuts junctions are 
not disordered. The legend provides the coordinates for each horizontal junction as in Fig. 1171 (a) p — 0.25, first shortcut 
realization. For 2.2 < a < 3.6, the horizontal junctions are weakly entrained, and for a > 3.6 the array is in the zero voltage 
state. The inset shows s v (a)/s v (0) versus a. The range of a values between the two vertical lines denote the region of weak 
entrainment. (b) p = 0.25, second shortcut realization. There is no evidence of entrainment up to a ~ 4.4, beyond which the 
array is in the zero voltage state. 

Setting the coefficients of the 0(e 1 ) terms in Eg. IA. 21 to zero gives: 

= 0c [dfai,j + 2d di7 0j ] + [doli.j + di 70j ] + S j - — R °-i+ s > ( A - 5 ) 

where we have divided by a factor of i c j. This result is written as 

#=0o7i,j = -2AA<9i7oj - <9i7oj - S ,j + — ^ R a .j+s- (A. 6) 

lc} 6=±1 
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FIG. 21: Cluster diagram for the time averaged-voltages, {vij) T , across the horizontal junctions, plotted versus coupling 
strength a for an array with M — 2 and N = 5. The bias current is is = 5, and the horizontal junctions are assigned critical 
currents randomly according to Eg. 1111 with A = 0.05. We assume periodic boundary conditions in the vertical direction. 
All shortcut junctions are identical to the uniform vertical junctions in the pristine array, i.e. the shortcuts junctions are 
not disordered. The legend provides the coordinates for each horizontal junction as in Fig. 1171 Shortcut junctions are added 
according to p = 1. 



Based on Eg. IA. 41 we can calculate the derivatives on the right side of Eg. IA.6I We find 
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dT dTi 



To 



/m,T 2) T 3 ) 



^-1=0 



because <pj is independent of To. Also 
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Next we look at the nonlinear terms on the right side of Eq. IA.6I 

'To , 



Sqj = sin 
Roj+s = sin 



lo. _ 



<t>3+5 ~ fij 



where in the limit of small disorder we have assumed To(l/i C j — l/i c .j + s) ~ 0. So Eo. IA.6l can be written as 



Pcdoli,j + 9 7ij = -<h<t>j - sin h <pj ) + M h 



(A.7) 



where Mj — sin [((f>j+8 — 4>j)/2] is a constant with respect to To- Temporarily ignoring the derivative di<f>j 

on the right side of Ea. IA.7I and using the trigonometric identity for the sine of a sum of two quantities, we find a 
particular solution to Ea. IA.7l of the form 



71J = M 3 To + Cf sin ( ?° ) + Dj 1 ' cos ( ^ 



(A. 
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where 
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(A.9) 
(A.10) 



The term MjTq in Eg. IA.8I represents a secular term that grows without bound as To — > oo. To remove this term 
from the solution we impose the condition 



Oo, + Mj = 0, 



which gives 



dcbj a v—* . 
T^r = — > sir 

dTx ^ 



S=±l 



i> j+s - 4>j 



(A.ll) 



This in turn gives a solution for 71 j that is Eg. I A. 81 without the secular term. Note that Eg. I A. Ill measures the rate 
of change of <pj, and hence 70. j, with respect to the slow time scale T\. 
Next, we look at all terms in Eg. IA. 21 that are 0(e 2 ): 

Pcdfa.j + d 72,j - -2/Wi7U - h (29 5 2 + df) 70., - 9i7i., - d 2l0 , 3 ~ S tJ + — ^ T i,j+s- (A. 12) 

tcj 6=±1 



Using the known results for jqj and jij to calculate the derivatives on the right side of Eq. IA.12I and using the 
expression for Sij and Ti^+s, means that Eo. IA.12l can be written (after some algebra) as 



where 



-d 2 <t>j + V 3 + Wj sin ( — J + Xj cos ( — J + Yj sin ( ?° ) + Z 3 cos ( — ) , (A.13) 
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where Cj and Z)j ' are given by Eos. IA.91 and lA.lOl As with the first order case, we want a solution to Eo. IA.131 
that does not have any secular terms. Therefore we must impose the condition 
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Based on Eg. IA. Ill it is possible to calculate d 2 <pj/dT 2 , which appears on the right side of Eg. IA.141 We find 
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where we approximated l/(i C j?cj±i) with 1/i 2 .,-- Note that V 2 0j = <j>j+i — 20j + (f>j—i. Substituting Eo. IA.151 into 
Eq. EH yields 
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(A.16) 



Our next step is to calculate the derivative djoj/dr with respect to the original dimensionless time variable 
t = f/is- We find 
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where use was made of the expressions e — l/iB and /3 C = is/3 c . It is convenient to define the quantity 
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(A.18) 



Physically, one can think of f2j as the angular frequency(average voltage) of oscillator(junction) j in the absence of 
coupling. Then Eq. IA.17I can be written as 



dr 



= % + r- E 



It is also useful to calculate the second derivative of 7oj 
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(A.19) 



where we made use of Eg. IA.151 It is common practice at this juncture to replace the symbol 70 j in Eo. IA.171 with 
(f>j , which we shall also do in Eq. IA. 191 

Finally, motivated by the structure of Eq. 1161 consider the combination of terms 

Substituting for the derivatives from Ea. lA.17l and lA.19l and dividing through by a factor of i C j results in the expression 
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which is Eq. in Sec. Ill Bl A straightforward but tedious continuation of the analysis to C(e 3 ) then leads to Eg. 1251 
in Sec. ITTBl 
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